An extended density matrix model applied to silicon-based terahertz quantum 

cascade lasers 
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Silicon-based terahertz quantum cascade lasers (QCLs) offer potential advantages over existing 
III-V devices. Although coherent electron transport effects are known to be important in QCLs, 
they have never been considered in Si-based device designs. We describe a density matrix transport 
model that is designed to be more general than those in previous studies and to require less a priori 
knowlege of electronic bandstructure, allowing its use in semi-automated design procedures. The 
basis of the model includes all states involved in interperiod transport, and our steady-state solution 
extends beyond the rotating-wave approximation by including DC and counter-propagating terms. 
We simulate the potential performance of bound-to-continuum Ge/SiGe QCLs and find that devices 
with 4-5-nm-thick barriers give the highest simulated optical gain. We also examine the effects of 
interdiffusion between Ge and SiGe layers; we show that if it is taken into account in the design, 
interdiffusion lengths of up to 1.5 nm do not significantly affect the simulated device performance. 

PACS numbers: 73.63.-b, 78.67.Pt, 05.60.Gg, 42.55.Px 
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I. INTRODUCTION 



Terahertz quantum cascade lasers (THz QCLs) are 
compact, coherent radiation sources in which electrons 
are transported through a periodic semiconductor het- 
erostructure, with a radiative transition in each period. 1 
Silicon-based THz QCLs may offer significant advan- 
tages over the existing III-V devices including the ab- 
sence of Reststrahlen absorption, which may allow emis- 
sion at frequencies above the current limit of 4.9 THz 
in GaAs/AlGaAs QCLs^ III-V devices require cryo- 
genic cooling (currently to < 200 K£) but the high ther- 
mal conductivity of Si and the absence of polar LO- 
phonon interactions could potentially overcome this lim- 
itation. Additionally, there is the prospect of leveraging 
the mature Si process technology, to create affordable 
integrated electrically-driven semiconductor THz lasers. 
Although mid-infrared^ and THz^ electroluminescence 
from p-type SiGe/Si quantum cascade structures has 
been achieved, lasing has not yet been demonstrated. 
In recent years, attention has switched to n-type de- 
vices owing to the simpler electron dispersion. We re- 
cently showed theoretically^ that the low effective mass 
and large usable conduction band offsets of the L-valleys 
in Ge/GeSi heterostructuresi^ could allow significantly 
higher gain, operating temperature and emission fre- 
quency range than equivalent Si/SiGe devices^ - — 

Several theoretical studies of Si-based QCLs 
exist A - -^— but none have accounted for coherent 
transport effects (i.e., quantum tunneling and inter- 
actions with optical fields). Although semiclassical 
scattering-transport models can give good agreement 
with experimental resultsp*^ they neglect tunneling 
across barriers, and can predict unrealistically large 
spikes in current density and gain when electrons scatter 
between spatially-extended subbandsJ^ By contrast, 
simplified density matrix (DM) models^— account 



for tunneling, in addition to scattering, and include 
the effect of the optical field on the electron dynamics. 
Additionally, DM models are much faster and less com- 
putationally demanding than full quantum mechanical 
simulations based on Green's functions^ - — To our best 
knowledge, all existing DM models of QCLs consider 
coherence between a reduced set of basis states including 
a single "injector" state (adjacent to a thick tunneling 
barrier) and a number of states in the next period of the 
QCL. This requires the manual selection of the injector 
state prior to simulating the device, and omits tunneling 
through the injection barrier from other states. This 
approach is not well-suited to semi-automated design 
proceduresj&S 8 - and can be problematic in bound-to- 
continuum (BTC) designs, in which multiple states 
may contribute to the tunneling current. Furthermore, 
these simplified models neglect coherences between the 
injector state and other states in the same period. In 
this work, we present a generalized DM model that 
reduces the requirement for a priori knowledge of the 
electronic bandstructure by including all states involved 
in interperiod transport, and by including contributions 
from both the optical field and the external DC bias in 
all density terms in our steady-state solution. 



Although the thickness of the injection barrier is of 
great importance in III-V QCls^ it has never been in- 
vestigated in Ge/GeSi devices. We therefore use our 
DM model in conjunction with a semi-automated de- 
sign algorithm 28 to investigate the influence of injec- 
tion barrier thickness upon the simulated device perfor- 
mance. Ge/GeSi interfaces also exhibit significantly more 
elemental interdiffusion than the GaAs/AlGaAs epitaxy 
used in existing QCLs. We use our DM model to as- 
sess the effect of interdiffusion on the simulated popula- 
tion inversion and gain and we compensate for the gain- 
reduction through design optimization. 
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II. THEORETICAL MODEL 
A. Bandstructure calculation 

The optically-active region of a QCL consists of a pe- 
riodic semiconductor heterostructure. In the DM model, 
the "injection barrier" that separates periods of the struc- 
ture is assumed to be sufficiently thick that interperiod 
transport is limited to quantum tunneling only. QCL pe- 
riods may be further subdivided into a number of mod- 
ules that are separated by thick tunneling barriers. 

We used a model-solid approximation to determine the 
conduction band profiles for the QCL structures in our 
simulations 3 ^ with the bandstructure parameters listed 
in Ref. |g. A self-consistent one-dimensional Poisson- 
Schrodinger solver was then used to locate the quasi- 
bound states within each module of the device, giving 
a total of TV subbands within each period. Localized 
wavefunctions tpi(z) with energies Ei were obtained, ac- 
cording to a 'tight-binding' scheme, by embedding each 
module of the structure between a pair of thick barriers, 
such that the amplitude of the evanescent waves decays 
to zero before reaching the edge of the simulation do- 
main. Here, the subscript i G [1 . . . N] denotes the index 
of each state in the period, in ascending order of energy. 
Intervalley mixing has previously been shown to yield 
negligibly small energy splitting for intcrsubband tran- 
sitions in Si-based heterostructures longer than a few 
nanometers ; 31 ! 32 and we therefore omit the effect from 
our model. 

As the QCL is a periodic heterostructure in an electric 
field, the states localized in other periods of the device 
are obtained by simple translations of these solutions in 
energy and space. The p th downstream period of the cas- 
cade therefore has states with energy Ef' =Ei- eFLp, 
where F is the applied electric field, e is the unit charge, 
and L is the length of a period. The corresponding wave- 
functions are given by ipl P \z) = ipi(z — Lp). The wave- 
function for an electron in the system is expressed in this 
basis as 



(i) 



( = 1 p=0 



where cf' is the weighting of each basis state and P is 
the number of periods contained in the model 
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B. Density matrix formulation 



coherences between these states and the injector describe 
interperiod tunneling transport. Some approaches sim- 
plify the calculation further by selecting only a subset of 
the states from the period after the barrier. The man- 
ual selection of an "injector" subband requires a priori 
knowledge of the electronic bandstructure, and is there- 
fore incompatible with semi-automated design optimiza- 
tion techniques. Even with a priori knowledge, the selec- 
tion of an injector state may not be obvious, particularly 
when the QCL is biased well away from subband align- 
ment; indeed, multiple channels for interperiod transport 
may exist. Also, this limited set of basis states does not 
include the injector subband in the second period. Al- 
though the effect on relaxation rates is included implicitly 
in the simulation, the calculation still does not account 
for any coherent interactions with this subband. Its role 
in intraperiod tunneling transport and its contribution 
to resonant optical transitions are therefore omitted. 

Here we describe a different DM model that uses all 
the subbands localized in three adjacent periods of the 
QCL as its basis. This method allows coherences to be 
calculated between all pairs of states in the central pe- 
riod, and allows interperiod tunneling (both in and out 
of the central period) to be determined without the need 
to select an injector subband manually. The resulting 
3N x 3./V density matrix may be expressed in block form 
as 



Pec Pcu Pcd\ 
P = | Puc Puu Pud 

yPDC PDU PDD / 



(2) 



where the subscripts U, C and D denote the upstream, 
center and downstream periods of the structure respec- 
tively. Each of these N x N blocks consists of density 
terms for pairs of states in the periods denoted by the 
block subscripts, for example, 



/pi,i ■•• Pi 



Pec 



N 



(3) 



\Pau ■ • • Pn, 



N J 



The density terms are unknown values to be solved, 
which represent the ensemble average of the weightings 
for the basis states pij = c*Cj. 

The Hamiltonian for the three-period system is written 
in block form as 



Hcc ^cu &CD 
^uc Huu ^UD 
,Qdc ^du Hdd, 



(4) 



Density matrix calculations rely on the selection of 
suitable basis states for coherent transport through the 
QCL. Conventional approaches use N basis states to 
model transport across a single injection barrier. A sin- 
gle "injector" subband is selected from the period up- 
stream of the barrier. The remaining N — 1 states are 
then selected from the period after the barrier, and the 



Here, the off-diagonal blocks contain the Rabi frequency 
terms for coupling between states in different periods of 
the QCL, which were calculated according to the scheme 
in Ref. 1 3 31 The diagonal blocks such as Hcc denote the 
Hamiltonian matrix for a single period of the structure. 
The elements of these single-period Hamiltonians are ei- 
ther the basis state energies (for diagonal elements), the 
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Rabi frequency terms (for pairs of states in different mod- 
ules) or the optical-coupling terms Hij = ZijA- mc for ra- 
diative transitions between states in response to incident 
light with the electric field A- mc = Ao exp(iwoi)- Here, the 
dipole matrix element terms are given by Zij — (ipi\z\i/jj). 

The time evolution of the density matrix is expressed 
by the Liouville equation 



dp 
dt 



dp 
di 



(5) 



relax 



where the last term is the matrix containing all relaxation 
and dephasing times. 

Although the Liouville equation for our three-period 
model contains a total of 9N 2 differential equations (com- 
pared with TV 2 for a single-period model), the transla- 
tional invariance of the system simplifies the calculation 
considerably. Firstly, the density terms may be trans- 
lated between blocks of the matrix, such that puu = 
Pec = Pdd, Pcd = Puc, and puc = Pcu- Similar 
translations may also be applied to the Rabi frequen- 
cies, dipole matrix elements, and relaxation times. The 
Liouville equation ([5]) can therefore be reduced to 3N 2 
independent differential equations, coming from the three 
blocks in the upper left corner of Eq. ((2} and (j4]). Sec- 
ondly, we apply the nearest-neighbor approximation, so 
that there is no coupling or scattering between states 
spaced by more than by one period. Therefore, the cou- 
pling constants between the first and third periods Qdu 
are set to zero in Eq. (QJ, and the same applies to the 
coherence and relaxation terms. These simplifications 
reduce the computational complexity of our model to a 
level comparable with the single-period model. 

The relaxation matrix can be written symbolically in 
the form 



[1/tcc, 1/t|i cc) (Vlcu) 
(1/t||uc) 



(6) 



where the dashes indicate that the corresponding N x N 
block is irrelevant owing to translational invariance. The 
relaxation terms in Eq. ([6]) determine both the linewidths 
of optical transitions and tunneling lifetimes. They con- 
tain contributions from the state lifetimes Tij, as well as 
'pure dephasing' times n ij for intraperiod transitions. 15 
The off-diagonal blocks in Eq. ([6]) (denoted T|| UC , etc) 
describe the pure dephasing for interperiod interactions. 

In this work, subband relaxation rates and tunneling 
dephasing rates were calculated by accounting for all rel- 
evant scattering processes^ acoustic and optical phonon 
scattering, intervalley scattering, ionized impurity and 
interface roughness scattering. The pure dephasing con- 
tributions are not straightforward to calculate, and a 
number of different schemes for their estimation have 
been proposed] 16 ' 20 ! 21 In the case of tunneling transitions, 
the prescription from Refs. and[2l| was employed. In 
case of optical transitions, however, we have simply set 
the linewidths to 2meV (as is typical for GaAs THz QCL 
structures) i 



C. Steady-state solution 

The harmonic balance method is a convenient ap- 
proach for determining a steady-state solution of the Li- 
ouville equation. Here, a simple steady-state functional 
form is assumed for each element of the density matrix 
depending on the states involved. Under the rotating- 
wave approximation (RWA), each element is assumed to 
contain only a single frequency harmonic. The diago- 
nal elements of the density matrix pij give the state 
populations, and are assumed to be constant-valued. 
The off-diagonal terms pij describe coherences between 
states. In the RWA approach, the steady-state forms 
of these terms are selected depending on whether they 
are optically- active or not. It is assumed that pairs of 
states within the same period with an energy separation 
Eij ~ ?ia;o±r will interact strongly with the driving opti- 
cal field (including a transition lincwidth F). The density 
matrix terms for optical transitions are assumed to be 
proportional to exp(iwoi) for Ei < Ej or cxp(— iwot) for 
Ei > Ej. All other off-diagonal elements, (i.e. for pairs of 
states with small energy separations, or in different mod- 
ules of the structure) are assumed to represent tunneling 
transport and are assumed to be constant valued. With 
the single dominant harmonic assumed for each p^j, the 
Liouville equation becomes a system of 37V 2 ordinary lin- 
ear equations. One of the equations is then replaced by 
the particle conservation law, Tr(p) = 1, and the system 
becomes inhomogeneous, with a unique solution. 

Although the RWA is useful for rapidly calculating the 
density matrix for a known system, it is incompatible 
with semi-automated design tools where the state ener- 
gies are not known a priori. The RWA also potentially 
omits multi-frequency effects on density matrix terms. 
In this work, therefore, we use an enhanced "non-RWA" 
solution method, that allows three harmonic terms to be 
included in the density matrix elements such that 

Pi,o = Pi,j exp(iw i) + P?j + Pi tj exp(-iwot), (7) 

where pfj, and pYf are unknown amplitudes for each 
of the harmonic components. As the Liouville equation 
is linear, each harmonic term may be treated indepen- 
dently, resulting in a system of up to 97V 2 linear equations 
(if every term contains all three harmonics), which still 
presents relatively low computational complexity. Simi- 
larly, in the light-field interaction terms in the Hamilto- 
nian both the exp(±iu;oi) components are retained. 

Previous studies indicate that the counter-propagating 
non-RWA terms may measurably affect the gain of lasers 
or the intensity-induced shift of resonance frequency in 
light-matter interactions, but are not very significant 
in near-resonant laser operationi2£~— Indeed, in all the 
cases considered in this work, we find that only one of 
the three harmonic amplitudes is significant in our simu- 
lations, which validates the use of the RWA when a priori 
knowledge of the state energies exists. 

The calculation can include thermal self-consistency 
(energy balance) by allowing the electron temperature T e 



4 



to be variable and requiring the total energy exchange be- 
tween the electron gas and the lattice to equal zeroJ^ In 
this work we did not include thermal balance within the 
density matrix equations, and have instead used the elec- 
tron temperature (assumed common to all subbands) as 
delivered by the rate equations model. 6 This is expected 
to be a reasonable approximation, since the tunneling 
contribution 3 ^ is not the major heat generating process 
in QCLs. 

D. Output parameters 

The current is calculated from the density matrix as 
j = Tr(,oJ), where 

J = ei[H,z]. (8) 

The current has both a time-independent (DC) compo- 
nent and a harmonic (AC) component that is induced by 
the optical field. The latter is used to find the complex 
permittivity e of the electron gas of the active medium, 
from 

_ d , . 

JAC = e ~fa A ™c (9) 

and then the gain from g w uj$ ■ Imag(e)/n r c, where n r is 
the refractive index and c is the speed of light in vacuum. 
Using very small values of A[ nc gives the small signal 
gain of the QCL active region (i.e. in the absence of gain 
saturation) . 

III. SIMULATION RESULTS 

We have recently reported the use of a semi-automated 
design optimization process to show that THz QCLs in 
the (001) Ge/GeSi material configuration yield substan- 
tially higher simulated gain than equivalent (111) and 
(001) Si/SiGe designs^ This method affords a system- 
atic comparison between different device design schemes 
and/or materials systems. In this work, we have used 
our optimization algorithm, in conjunction with the DM 
model described above, to identify viable designs for 
bound-to-continuum (BTC) Ge/SiGe QCLs, while ac- 
counting for coherent transport effects. Previous DM 
models of III-V BTC QCLs have reproduced experimen- 
tal results by including only the upper laser level (ULL) 
and a subset of miniband states in their basis set, of 
which one is designated as the injector^ However, in 
practice multiple subbands may contribute significantly 
to interperiod tunneling in BTC QCLs. Our extended 
non-RWA DM model avoids the need for a priori selec- 
tion of the optically-coupled transition and the injector 
subband, allowing us to use the semi-automated design 
approach described in Ref. l28l . Furthermore, our model 
explicitly includes all possible interperiod tunneling path- 
ways. 




Position [nm] 



FIG. 1. Two periods of the bandstructure of the 
six- well BTC QCL design which serves as a template 
for our investigations. The thicknesses of the lay- 
ers are as follows (starting from the injection bar- 
rier): 4.3/5.0/1.2/14.3/1.3/12.9/1.6 /7.9/1.8/7.1/2.5/7.8 
where bold text denotes the Sio.15Geo.85 barriers and regu- 
lar text denotes Ge wells. The underlined layers are doped to 
provide a sheet doping density of 8 x 10 10 cm -2 . The upper 
and lower laser level in each period of the structure are shown 
as bold, and dashed-bold lines respectively. 



Figure [T] shows two periods of the bandstructure of a 
4 THz six- well BTC QCL design which serves as a tem- 
plate for our investigation of device performance. For our 
DM simulations, each period of the active-region struc- 
ture was modeled as a single module. All the devices 
simulated in this work were derived from this template 
by systematically adjusting a single element of the struc- 
ture (i.e. either the injection barrier thickness or the in- 
tcrdiffusion length) . In all lattice temperature of 

4K was used. The electron temperature was fixed at a 
value of 100 K, which was obtained from a semiclassical 
simulation of the design template. 

As Si and Ge are not lattice matched, the proposed 
QCL structures must be grown on a relaxed buffer or 
"virtual substrate" . Such substrates can be achieved ei- 
ther by growth of a linearly graded alloy layer, from pure 
silicon up to the desired buffer composition ) 3 ^ 39 or by 
growth of a Ge seed layer directly on a silicon wafer, fol- 
lowed by reverse linear grading from pure Ge down to the 
buffer composition i 40 i 41 The relaxed buffer composition 
is calculated to achieve "strain symmetrization" through- 
out the QCL structure^ whereby the compressive stress 
in the Ge wells is balanced by the tensile stress in the bar- 
riers, yielding zero net stress over each period. For all the 
cases presented below, the optimum buffer composition 
was found to be Geo.97Sio.03- 
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FIG. 2. (a) Simulated current density and (b) gain spectra of optimized BTC QCLs for different injection barrier thicknesses. 



An indicative figure for the lasing threshold is included for reference, based on calculations of the waveguide losses in Ref. [43 



A. Injection barrier thickness 

The thickness of the injection barrier is known to be 
an important parameter in III-V THz QCLs as it can 
significantly affect the performance of devices^ If the 
injection barrier is too thin, selectivity of injection into 
the upper laser level is poor; if it is too thick, the tunnel- 
ing rate through the barrier is small and efficient injection 
cannot be achieved. Semiclassical rate-equation models 
of charge injection in THz QCLs lack sensitivity to the 
injection barrier thickness, and a model that accounts for 
coherent effects is therefore required^ Here, we system- 
atically vary the thickness of the injection barrier in the 
design template, and use our DM simulation to deter- 
mine the gain and current density. We use the genetic 
algorithm described in Ref. [28| to maximize the simulated 
gain of the device at 4 THz by varying the thickness of 
each of the other layers in the structure, and the applied 
electric held. 

In all the cases considered here, the optimized layer 
widths (excluding the injection barrier) were found to 
be identical (to angstrom precision) to those of the tem- 
plate. This can be understood by recalling that the de- 
coupled wavefunctions in the DM model are found by em- 
bedding the active-region module between thick barriers. 
Therefore, the layer-width optimization procedure only 
directly affects the single-period Hamiltonian, and has a 
very much weaker influence on the interperiod coupling 
(by slightly adjusting the Rabi frequencies). As such, the 
results below can be seen to solely represent the effect of 
the injection barrier thickness upon the device perfor- 
mance without including any contribution from changes 
in the bandstructure within the period. 

Figure [5a| shows the simulated current density as a 
function of applied electric field for optimized devices 
with 2, 3, 4, 5, and 6-nm-thick injection barriers. Interpe- 
riod scattering between spatially-extended wavefunctions 



is avoided in the DM model, and the simulated current 
density is therefore a smoothly-varying function of bias 
in all cases. The alignment bias for the device is ap- 
proximately 4kV/cm for all five structures, and we see 
that the current density at this bias decreases monoton- 
ically as the thickness of the injection barrier increases. 
The gain spectrum for each structure at the alignment 
bias is shown in Fig. I2bl The variation in the magnitude 
of the peak gain with injection barrier thickness is not 
monotonic, and there is an optimum thickness at around 
4-5 nm, at which a gain of around 90 cm^ 1 is predicted. 
The difference in gain spectrum between the device with 
2-nm-thick injection barriers and the other structures is 
caused principally by the reduction of injection selectiv- 
ity into the ULL. It is important to note, however, that 
the "injection barrier" in this structure is thinner than 
the 2.5-nm-thick barrier at the end of the module, and 
the chosen subdivision of modules for tunneling transport 
in the DM model is likely to be unrealistic. 



B. Interdiffusion compensation 

Epitaxial Ge/SiGe heterostructures have been re- 
ported to show significant interdiffusion between the pure 
and alloy semiconductor layers, with typical characteris- 
tic interdiffusion lengths estimated to be of the order of 
1-2 nmj^fi This leads to significant changes in the band- 
structure and scattering lifetimes , 46 ' 47 and can therefore 
degrade the performance of devices. In this section, we 
investigate the impact of interdiffusion on the QCL gain 
and we attempt to recover the lost performance through 
design optimization. 

We account for the effects of interdiffusion by apply- 
ing a Gaussian annealing model to the alloy composition 
profile— such that the alloy fraction across the interface 
is described by a Gauss error function with the interdif- 
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fusion length Ld as a size parameter. Figure |3] shows the 
bandstructure of the QCL design template with 1 and 
2nm interdiffusion included. The interdiffusion causes 
the barriers to become reduced in height and the shape 
of the quantum wells becomes distorted, with the tops 
being widened and the bottoms being narrowed. For in- 
terdiffusion lengths of 2 nm, the thinner barriers near the 
optically active wells are significantly reduced in height 
and the ULL is poorly confined. 

Fig. I4al shows the simulated gain spectra for structures 
with Ld in the range 0-2 nm at the respective operat- 
ing bias for each device. Here, the interdiffusion has 
been included without subsequently optimizing the de- 
sign. There is no significant drop in the peak gain when 
1 nm interdiffusion length is included. However, for larger 
interdiffusion lengths the simulated peak gain is reduced 
considerably, and there is no simulated gain at all for the 
structure with Ld = 2 nm. 

We applied the design optimization algorithm to each 
of the diffuse QCL structures, and the resulting gain spec- 
tra are shown in Fig. I4bl It can be seen that the gain 
has been fully recovered for structures with interdiffusion 
lengths up to 1.5 nm. This highlights the importance 
of being able to characterize the interdiffusion length in 
these system: so long as this is known, and so long as it 
is 1.5 nm or less, our results indicate that it can be taken 
into account in the design process. 

For interdiffusion lengths of 2 nm or more, the gain 
cannot be recovered. There are two reasons for this: 
first, the ULL is no longer confined by the thin barriers 
in the structure, which leads to a large spatial overlap 
with the miniband states, and hence a loss of population 
inversion. Second, the interdiffusion introduces Si into 
the nominally pure Ge well regions, leading to a large 
increase in alloy disorder scattering rates4£ As such, ad- 
ditional rapid scattering pathways are introduced to the 
system, leading to rapid depopulation of the ULL. By 



way of comparison, we have previously shown that the 
total scattering rate within a Si/Ge/Si quantum well in- 
creases by 50% at an interdiffusion length of 1.21 nm4£ 



IV. CONCLUSION 

We have investigated coherent transport effects in a Si- 
based THz QCL through the use of an extended density 
matrix model that includes in its basis all subbands that 
are involved in interperiod transport. Our use of the 
non-rotating wave approximation allows a steady-state 
solution to the Liouville equation without a priori knowl- 
edge of the bandstructure of the device. In all cases, the 
non-RWA solution yielded a single strongly-dominant fre- 
quency component in each density term, indicating that 
it would be in good agreement with an equivalent RWA 
model. Although we have used our generalised non-RWA 
model to analyze coherent effects in Si-based QCLs, it is 
equally applicable to III-V QCL structures. 

We have coupled our model with a semi-automated 
QCL design algorithm, and have shown that the optimum 
injection barrier thickness for Si-based BTC THz QCL 
structures is in the range 4-5 nm, and we predict peak 
gain values of ~ 90 cm -1 at a lattice temperature of 4K. 
We have also studied the effect of interdiffusion between 
the Ge and GeSi layers, and found that it is possible 
to compensate for interdiffusion effects through design 
optimization up to a limit of Ld ~ 1.5 nm. 
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